\documentclass{ctexart}
\usepackage{listings}%插入代码
\usepackage{geometry}%设置页面大小边距等
\usepackage{graphicx}%插入图片
\usepackage{amssymb}%为了用\mathbb
\usepackage{amsmath}%数学方程的显示
\usepackage{listings}%插入代码
\usepackage{fancyhdr}%设置页眉页脚
\usepackage{lastpage}%总页数
\usepackage{hyperref}%引用网页
\usepackage{xcolor}
\usepackage{tikz}
\usepackage{float}


\geometry{a4paper,left=2cm,right=2cm,top=2cm,bottom=2cm}%一定要放在前面！
\pagestyle{fancy}%设置页眉页脚
\lhead{陈冠宇\ 3200102033}%页眉左
\chead{Numerical Methods for Differential Equations}%页眉中
\rhead{Project2}%章节信息
\cfoot{\thepage/\pageref{LastPage}}%当前页，记得调用前文提到的宏包
\lfoot{Zhejiang University}
\rfoot{School of Mathematical Sciences}
\renewcommand{\headrulewidth}{0.1mm}%页眉线宽，设为0可以去页眉线
\renewcommand{\footrulewidth}{0.1mm}%页脚线宽，设为0可以去页眉线
\setlength{\headwidth}{\textwidth}

\hypersetup{%设置网页链接颜色等
    colorlinks=true,%链接将会有颜色，默认是红色
    linkcolor=blue,%内部链接，那些由交叉引用生成的链接将会变为蓝色（blue）
    filecolor=magenta,%链接到本地文件的链接将会变为洋红色（magenta）
    urlcolor=blue,%链接到网站的链接将会变为蓝绿色（cyan）
    }

\newtheorem{theorem}{Theorem}
\newtheorem{proof}{Proof}
\newtheorem{solution}{Solution:}

\title{\textbf{微分方程数值解第二次大作业}}
\date{\today}

\begin{document}
\section*{Project2}
测试方法project2文件夹中make即可，之后的可执行文件在bin文件夹中。
\subsection*{一维区域中的BVP(Assignment.A)}
由于一维条件的特殊性，报告中只展示Mixed边界条件。最大迭代次数为100, 最小误差为$\epsilon = 10^{-8}, \nu_1 = \nu_2 = 10$.
\subsection*{\textbf{$f(x) = e^{sin(x)}$}}
Combination 1: Restriction: injection; Interpolation: Linear;

当$n = 32$时
\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &4.96132e-05 &离散解误差： &1.36937e-05\\ \hline
    V-cycle: 2th: &解析解误差： &4.94151e-05 &离散解误差： &4.5258e-10\\ \hline
    FMG: 1th:&解析解误差： &4.87519e-05 &离散解误差： &7.40048e-07\\ \hline
    \end{tabular}
\end{table}

Similarly, n = 64, 128, 256. 由于篇幅有限，再将n=256作展示：

\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &1.373e-05 &离散解误差： &1.37034e-05\\ \hline
    V-cycle: 2th: &解析解误差： &7.73163e-07 &离散解误差： &6.32893e-10 \\ \hline
    FMG: 1th:&解析解误差： &3.2191e-07 &离散解误差： &8.93316e-07\\ \hline
    \end{tabular}
\end{table}

Combination 2: Restriction: full weighting; Interpolation: quadratic;

由于测试过程相似，报告中仅展示$n = 256$时，

\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.0641714 &离散解误差： &0.0641713\\ \hline
    V-cycle: 2th: &解析解误差： &0.00148234 &离散解误差： &0.00148231\\ \hline
    V-cycle: 3th: &解析解误差： &2.95913e-05 &离散解误差： &2.95647e-05\\ \hline
    V-cycle: 4th: &解析解误差： &7.74182e-07 &离散解误差： &5.07988e-07\\ \hline
    V-cycle: 5th: &解析解误差： &7.72717e-07 &离散解误差： &7.45182e-09\\ \hline
    FMG: 1th: &解析解误差： &0.000239103 &离散解误差： &0.000239085 \\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/11.png}
  \caption{$f(x) = e^{sin(x)}$}
  \label{fig:example}
\end{figure}

\subsection*{\textbf{$f(x) = x^2 - x$}}
同样的测试方式，我们仅将$n = 256$时作展示：(Combination 同上)

\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.00392375 &离散解误差： &0.00392375\\ \hline
    V-cycle: 2th: &解析解误差： &8.39099e-05 &离散解误差： &8.39099e-05\\ \hline
    V-cycle: 3th: &解析解误差： &2.10351e-06 &离散解误差： &2.10351e-06\\ \hline
    V-cycle: 4th: &解析解误差： &4.39843e-08 &离散解误差： &4.39843e-08\\ \hline
    V-cycle: 5th: &解析解误差： &7.49741e-10 &离散解误差： &7.49741e-10\\ \hline
    FMG: 1th:&解析解误差： &1.61356e-06 &离散解误差： &1.61356e-06 \\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/12.png}
  \caption{$f(x) = x^2 - x$}
  \label{fig:example}
\end{figure}


\subsection*{\textbf{$f(x) =\frac{1}{x+1}$}}
同样的测试方式，我们仅将$n = 256$时作展示：(Combination 同上)

\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.0149794 &离散解误差： &0.0149794\\ \hline
    V-cycle: 2th: &解析解误差： &0.00038608 &离散解误差： &0.000386108\\ \hline
    V-cycle: 3th:& 解析解误差： &8.53876e-06 &离散解误差： &8.56596e-06\\ \hline
    V-cycle: 4th: &解析解误差： &7.30116e-07 &离散解误差： &1.62954e-07\\ \hline
    V-cycle: 5th: &解析解误差： &7.3102e-07 &离散解误差： &2.6375e-09\\ \hline
    FMG: 1th: &解析解误差： &5.23781e-05 &离散解误差： &5.23917e-05 \\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/13.png}
  \caption{$f(x) =\frac{1}{x+1}$}
  \label{fig:example}
\end{figure}


\subsection*{二维区域中的BVP(Assignment.B)}
二维中选择Dirichlet边界进行测试，最大迭代次数为100, 最小误差为$\epsilon = 10^{-8}, \nu_1 = \nu_2 = 10$.
\subsection*{\textbf{$f(x) = e^{sin(\pi * x)}$}}
Combination 1: Restriction: injection; Interpolation: Linear;
由于在设计中使用了稠密矩阵，当网格大小很大时，会出现数组越界的情况，当n = 32时：
\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.0813769 &离散解误差： &0.0813534 \\ \hline
	V-cycle: 2th: &解析解误差： &0.00170508 &离散解误差： &0.00167242\\ \hline
	V-cycle: 3th: &解析解误差： &6.72193e-05 &离散解误差： &3.35533e-05\\ \hline
	V-cycle: 4th: &解析解误差： &3.50564e-05 &离散解误差： &6.54272e-07\\ \hline
	V-cycle: 5th: &解析解误差： &3.44618e-05 &离散解误差： &1.26363e-08\\ \hline
	V-cycle: 6th: &解析解误差： &3.44505e-05 &离散解误差： &2.43026e-10\\ \hline
    FMG: 1th:&~&~&离散解误差：&2.73519e-05\\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/21.png}
  \caption{$f(x) = e^{sin(\pi * x)}$}
  \label{fig:example}
\end{figure}


\subsection*{\textbf{$f(x) = x^2 + y^2$}}
Combination 1: Restriction: injection; Interpolation: Linear;
由于在设计中使用了稠密矩阵，当网格大小很大时，会出现数组越界的情况，当n = 32时：
\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.0218465 &离散解误差： &0.0218465\\ \hline
	V-cycle: 2th: &解析解误差： &0.000379184 &离散解误差： &0.000379184\\ \hline
	V-cycle: 3th: &解析解误差： &7.51058e-06 &离散解误差： &7.51058e-06\\ \hline
	V-cycle: 4th: &解析解误差： &1.45623e-07 &离散解误差： &1.45623e-07\\ \hline
	V-cycle: 5th: &解析解误差： &2.80741e-09 &离散解误差： &2.80748e-09\\ \hline
    FMG: 1th:&~&~&离散解误差：&3.36207e-05\\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/22.png}
  \caption{$f(x) = x^2 + y^2$}
  \label{fig:example}
\end{figure}

\subsection*{\textbf{$f(x) = sin(x+y)$}}
Combination 1: Restriction: injection; Interpolation: Linear;
由于在设计中使用了稠密矩阵，当网格大小很大时，会出现数组越界的情况，当n = 32时：
\begin{table}[!ht]
    \centering
    \begin{tabular}{|l|l|l|l|l|}
    \hline
    V-cycle: 1th: &解析解误差： &0.0225398 &离散解误差： &0.0225496\\ \hline
	V-cycle: 2th: &解析解误差： &0.000471652 &离散解误差： &0.000481488\\ \hline
	V-cycle: 3th:& 解析解误差： &1.85838e-06& 离散解误差： &9.6433e-06\\ \hline
	V-cycle: 4th: &解析解误差： &9.64907e-06 &离散解误差： &1.88533e-07\\ \hline
	V-cycle: 5th: &解析解误差： &9.8326e-06 &离散解误差： &3.64257e-09\\ \hline
	FMG: 1th:&~&~&离散解误差：1.32881e-05\\ \hline
    \end{tabular}
\end{table}

图像如下：
\begin{figure}[H]
  \centering
  \includegraphics[width=0.5\textwidth]{pic/23.png}
  \caption{$f(x) = sin(x+y)$}
  \label{fig:example}
\end{figure}

\section*{计算时间和收敛阶数对比}
由于使用稠密矩阵缘故，算的很慢。

\section*{理论分析}
\subsection*{Question}
\par 考虑二维Dirichlet边值条件的BVP：
\begin{equation}
	\begin{cases}
		-\Delta u = f, & in\ \Omega:=(0,1)^2 \\
		u = 0, & on\ \partial\Omega
	\end{cases}
	\label{BVP}
\end{equation}

	
\subsection*{Relaxation}
\subsection*{Thm}
	二维BVP求解中的线性方程组系数矩阵$A_{2D} = I\otimes A + A\otimes I$，其中$A$是一维Dirichlet边值条件BVP线性系统的系数矩阵.并且$A_{2D}$的特征值和特征向量是
	$$\lambda_{ij}=\lambda_i+\lambda_j,\quad \mathbf{W}_{ij}=\mathbf{w}_j\otimes\mathbf{w}_i$$
	其中$\lambda_i$和$\lambda_j$是A的特征值，$\mathbf{w}_i$和$\mathbf{w}_j$是对应的特征向量.
\begin{proof}
	这直接由讲义中的Theorem 7.48和Theorem 7.51得到.
\end{proof}

\subsection*{Thm}
	\eqref{BVP}对应的线性系统的加权Jacobi迭代矩阵是
	$$T_{\omega} = (1-\omega)I+\omega D^{-1}(L+U) = I-\frac{\omega h^2}{4}A_{2D}$$
	其特征向量就是$A_{2D}$的特征向量，对应的特征值是
	$$\lambda_{ij} = 1-\omega\sin^2(\frac{i\pi}{2n})-\omega\sin^2(\frac{j\pi}{2n})$$

\begin{proof}
	这由讲义Lemma 9.16直接得到.
\end{proof}

\subsection*{Restriction and prolongation}
\subsection*{Lemma}
	二维网格的full weighting算子和quadratic算子满足
	$$I_{h(2D)}^{2h} = I^{2h}_{h}\otimes I^{2h}_{h},\quad I_{2h(2D)}^{h}=I^{h}_{2h}\otimes I^{h}_{2h}$$
	其中$I_{h}^{2h}$和$I_{2h}^{h}$分别时一维的full weighting算子和quadratic算子.

    \subsection*{Example}
	考虑在$n=3$的情况，我们有二维full weighting算子为
	\begin{align*}
		I^{2h}_{h(2D)} &= \frac{1}{16}\begin{bmatrix}
			1 & 1 & -1 & 1 & 8 & 1 &  1 & 1 &  1 
		\end{bmatrix} 
	\end{align*}


\subsection*{Two-grid correction}
\subsection*{Lemma}
	二维的Two-grid correction作用在误差向量上的迭代矩阵是
	$$TG=T_{\omega}^{\nu_2}[I_{2D}-I^{h}_{2h(2D)}(A^{2h}_{2D})^{-1}I^{2h}_{h(2D)}A^h_{2D}]T_{\omega}^{\nu_1}$$

\begin{proof}
	这由讲义Lemma 9.31推广得到.
\end{proof}

\subsection*{The spectral picture}
\subsection*{Lemma}
	full weighting restriction算子作用在$A_{2D}^{h}$的特征向量$\mathbf{W}_{ij}^{h}$上得到
	\begin{align*}
		I^{2h}_{h(2D)}\mathbf{W}_{ij}^{h} := c_ic_j\mathbf{W}_{ij}^{2h} = \cos^2\frac{i\pi}{2n}\cos^2\frac{j\pi}{2n}\mathbf{W}_{ij}^{2h}
	\end{align*}

\begin{proof}
	\begin{align*}
		I^{2h}_{h(2D)}\mathbf{W}_{ij}^{h} &= (I^{2h}_h\otimes I^{2h}_h)\cdot(\mathbf{w}_j^h\otimes\mathbf{w}_i^h) \\
		&= (I^{2h}_{h}\mathbf{w}_{j}^{h})\otimes(I^{2h}_{h}\mathbf{w}_{i}^{h}) \\
		&= (c_i\mathbf{w}_{j}^{2h})\otimes(c_j\mathbf{w}_{i}^{2h}) \\
		&= c_ic_j\mathbf{W}_{ij}^{2h}
	\end{align*}
\end{proof}

\subsection*{Lemma}
	linear interpolation算子作用在$A_{2D}^{h}$的特征向量$\mathbf{W}_{ij}^{h}$上得到
	$$I^h_{2h(2D)}\mathbf{W}_{ij}^{2h} = c_ic_j\mathbf{W}_{ij}^h-c_is_j\mathbf{W}_{ij'}^h-c_js_i\mathbf{W}_{i'j}^h+s_is_j\mathbf{W}_{i'j'}^h$$
	其中，$c_i = \cos^2\frac{i\pi}{2n}$，$s_i=\sin^2\frac{i\pi}{2n}$，$i'=n-i$，$j'=n-j$.

\begin{proof}
	\begin{align*}
		I^h_{2h(2D)}\mathbf{W}_{ij}^{2h} &= (I^h_{2h}\otimes I^h_{2h})\cdot(\mathbf{w}_{j}^{2h}\otimes\mathbf{w}_{i}^{2h}) \\
		&= (I^h_{2h}\mathbf{w}_{j}^{2h})\otimes(I^h_{2h}\mathbf{w}_{i}^{2h}) \\
		&= (c_j\mathbf{w}_{j}^{2h}-s_j\mathbf{w}_{j'}^h)\otimes(c_i\mathbf{w}_{i}^{2h}-s_i\mathbf{w}_{i'}^h) \\
		&= c_ic_j\mathbf{W}_{ij}^h-c_is_j\mathbf{W}_{ij'}^h-c_js_i\mathbf{W}_{i'j}^h+s_is_j\mathbf{W}_{i'j'}^h
	\end{align*}
\end{proof}

\subsection*{Thm}
	Two-grid correction作用在子空间$span\{\mathbf{W}_{ij}^{h},\mathbf{W}_{i'j}^{h},\mathbf{W}_{ij'}^{h},\mathbf{W}_{i'j'}^{h}\}$上是不变的，且
	\begin{align*}
		TG\mathbf{W}_{ij}^{h} &= \lambda_{ij}^{\nu_1+\nu_2}\left(1-\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_ic_j\right)\mathbf{W}_{ij} + \lambda_{ij}^{\nu_1}\lambda_{ij'}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_is_j\mathbf{W}_{ij'}^{h} \\ &+ \lambda_{ij}^{\nu_1}\lambda_{i'j}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_js_i\mathbf{W}_{i'j}^{h} - \lambda_{ij}^{\nu_1}\lambda_{i'j'}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}s_is_j\mathbf{W}_{i'j'}^{h}
	\end{align*}
	其中，$\lambda_{ij}=1-\omega\sin^2(\frac{i\pi}{2n})-\omega\sin^2(\frac{j\pi}{2n})$是$T_{\omega}$的特征值，$\chi_i=\sin^2\frac{i\pi}{n}$，$\chi_j=\sin^2\frac{j\pi}{n}$.

\begin{proof}
	首先不考虑relaxation过程，即假设$\nu_1=\nu_2=0$，
	\begin{align*}
		&A^h_{2D}\mathbf{W}_{ij}^h=\frac{4}{h^2}(s_i+s_j)\mathbf{W}_{ij}^h \\
		\implies& I^{2h}_{h(2D)}A^h_{2D}\mathbf{W}_{ij}^h=\frac{4}{h^2}(s_i+s_j)c_ic_j\mathbf{W}_{ij}^{2h} \\
		\implies& (A^{2h}_{2D})^{-1}I^{2h}_{h(2D)}A^h_{2D}\mathbf{W}_{ij}^h=\frac{h^2}{\chi_i+\chi_j}\frac{4}{h^2}(s_i+s_j)c_ic_j\mathbf{W}_{ij}^{2h}=\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}\mathbf{W}_{ij}^{2h} \\
		\implies& I_{2h(2D)}^{h}(A^{2h}_{2D})^{-1}I^{2h}_{h(2D)}A^h_{2D}\mathbf{W}_{ij}^h=\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}[c_ic_j\mathbf{W}_{ij}^h-c_is_j\mathbf{W}_{ij'}^h-c_js_i\mathbf{W}_{i'j}^h+s_is_j\mathbf{W}_{i'j'}^h] \\
		\implies& [I-I_{2h(2D)}^{h}(A^{2h}_{2D})^{-1}I^{2h}_{h(2D)}A^h_{2D}]\mathbf{W}_{ij}^h=\left(1-\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_ic_j\right)\mathbf{W}_{ij}^h \\
		&+ \frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_is_j\mathbf{W}_{ij'}^{h} + \frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_js_i\mathbf{W}_{i'j}^{h} - \frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}s_is_j\mathbf{W}_{i'j'}^{h}
	\end{align*}
	其中第二步和第四步来自上述两个引理，考虑上smoothing的过程得到：
	\begin{align*}
	TG\mathbf{W}_{ij}^{h} &= \lambda_{ij}^{\nu_1+\nu_2}\left(1-\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_ic_j\right)\mathbf{W}_{ij}^h + \lambda_{ij}^{\nu_1}\lambda_{ij'}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_is_j\mathbf{W}_{ij'}^{h} \\ &+ \lambda_{ij}^{\nu_1}\lambda_{i'j}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}c_js_i\mathbf{W}_{i'j}^{h} - \lambda_{ij}^{\nu_1}\lambda_{i'j'}^{\nu_2}\frac{\chi_ic_j+\chi_jc_i}{\chi_i+\chi_j}s_is_j\mathbf{W}_{i'j'}^{h}
	\end{align*}
\end{proof}

\subsection*{Thm}
	类似上述讨论，$TG$作用在$\mathbf{W}_{ij'}^h$，$\mathbf{W}_{i'j}^h$，$\mathbf{W}_{i'j'}^h$上也有四项系数充分小于$1$，即：
	\begin{align*}
	TG\begin{bmatrix}
		\mathbf{W}_{ij}^h \\
		\mathbf{W}_{ij'}^h \\
		\mathbf{W}_{i'j}^h \\
		\mathbf{W}_{i'j'}^h
	\end{bmatrix}=\begin{bmatrix}
		k_{11} & k_{12} & k_{13} & k_{14} \\
		k_{21} & k_{22} & k_{23} & k_{24} \\
		k_{31} & k_{32} & k_{33} & k_{34} \\
		k_{41} & k_{42} & k_{43} & k_{44} 
	\end{bmatrix} \begin{bmatrix}
		\mathbf{W}_{ij}^h \\
		\mathbf{W}_{ij'}^h \\
		\mathbf{W}_{i'j}^h \\
		\mathbf{W}_{i'j'}^h
	\end{bmatrix}
	\end{align*}
	其中，$k_{ij}$都是充分小于$1$的数.

综上，说明了二维多重网格方法的收敛性.

\end{document} 